% This code simulates the model for EIS ~= 1 under RE. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function SIMDATA = Sim_EIS_RE(l_q , l_q_burn , N_sim , psi , c_0 , Var_Names)

% INPUT:
% -- l_q: Number of sample periods
% -- l_q_burn: Number of burn-in periods
% -- N_sim: Number of simulations
% -- psi: EIS parameter
% -- c_0: starting value of log real per capita consumption level
% -- Var_Names: names of variables to be stored
%
% OUTPUT:
% -- SIMDATA: stored simulations as a cell

%% 1. LOAD PARAMETERS
[I , J] = param_EIS_RE(psi);

%% 2. SIMULATIONS
% additional parameters for simulation
T = l_q;                %Number of sample periods
T_burn = l_q_burn;            %Number of burn-in periods
T_total = T_burn + T;     %Total periods including burn-in periods for simulation
d_0 = c_0 + I.d_c_mean;   %Starting value of log dividend

%% 2.1 Simulate the economy
% Simulate quarterly iid consumption growth, column for each simulation
cg_full = normrnd(I.mu , I.sigma , T_total , N_sim); 

% Simulate quarterly shocks to dividend growth
dg_shock = normrnd(0 , I.sigma_d , T_total , N_sim); 

% Log real consumption level from 1 to T_total
c_full = c_0 + cumsum(cg_full); 

% Log real dividend level from 1 to T_total
d_full = zeros(T_total , N_sim); 

d_full(1 , :) = I.lambda * cg_full(1 , :) + (1 - I.alpha) * d_0 + I.alpha * c_0 + I.alpha * I.mu_dc + dg_shock(1 , :);

for k = 2 : T_total
  d_full(k , :) = I.lambda * cg_full(k , :) + (1 - I.alpha) * d_full(k - 1 , :) + I.alpha * c_full(k - 1 , :) + I.alpha * I.mu_dc + dg_shock(k , :);
end

% Log dividend growth from 1 to T_total
dg_full = [d_full(1 , :) - d_0 ; d_full(2 : end , :) - d_full(1 : end - 1 , :)]; 

% Four-quarter change
cg_yoy_full = (c_full(5 : end , :) - c_full(1 : end - 4 , :))/4;
dg_yoy_full = (d_full(5 : end , :) - d_full(1 : end - 4 , :))/4;

% Construct State Variable \tilde{\mu}
tmu_full = I.mu * ones(T_total , N_sim);

% Calculate risk-free rate
rf_full = -I.tmu_m_s + tmu_full/I.psi - 0.5 * I.gamma^2 * I.sigma^2;

%% 2.2 Calculate price-dividend ratio and price
pd_full = J.Ad_0 + J.Ad_2 * (d_full - c_full);

P_full = exp(pd_full + d_full);

% Log return implied by log-linearization
ret_ll_full = J.kappad_0 ...
            + J.kappad_1 * pd_full(2 : end , :) - pd_full(1 : end - 1 , :) ...
            + dg_full(2 : end , :);
        
% Log return calculated by price and dividend
ret_full = log( P_full(2 : end , :) + exp(d_full(2 : end , :)) ) ... % log(P_{t+1} + D_{t+1})
         - log( P_full(1 : end - 1 , :) ); % - log(P_t)
     
% Excess returns     
EXRET_full = exp(ret_full) - exp(rf_full(1 : end - 1 , :)); % Net excess returns
EXRET_ll_full = exp(ret_ll_full) - exp(rf_full(1 : end - 1 , :)); % Net excess returns

exret_full = ret_full - rf_full(1:end-1,:); %log excess returns
exret_ll_full = ret_ll_full - rf_full(1:end-1,:); %log excess returns

% Trailing 4-quarter sum of dividends
D_ttm_full = movsum(exp(d_full) , 4 , 'Endpoints' , 'discard');

PD_ttm_full = P_full(4 : end , :) ./ D_ttm_full; %For the last t periods
pd_ttm_full = log(PD_ttm_full);

%% 3. SIMULATED DATA

%% 3.1 Drop burn-in periods
% All data are re-labled from 1 to T as realizations. 
% Be careful with risk-free rate

%-------------------------------------
% Endowment Process
%-------------------------------------
cg = cg_full(end - T + 1 : end , :);
dg = dg_full(end - T + 1 : end , :);
c  =  c_full(end - T + 1 : end , :);
d  =  d_full(end - T + 1 : end , :);
cg_yoy = cg_yoy_full(end - T + 1 : end , :);
dg_yoy = dg_yoy_full(end - T + 1 : end , :);

%-------------------------------------
% State variable proxies
%-------------------------------------
%tmu  =  tmu_full(end - T + 1 : end , :);
tmu = nan(T , N_sim);
tmud = nan(T , N_sim);
tmur = nan(T , N_sim);
LTG = nan(T , N_sim); 
Y1G = nan(T , N_sim);

%-------------------------------------
% Price Quantity
%-------------------------------------
rf = rf_full(end - T + 1 : end , :);

P = P_full(end - T + 1 : end , :);
D_ttm = D_ttm_full(end - T + 1 : end , :);
pd_ttm = pd_ttm_full(end - T + 1 : end , :);
pd = pd_full(end - T + 1 : end , :);

ret = ret_full(end - T + 1 : end , :);
ret_ll = ret_ll_full(end - T + 1 : end , :);

EXRET = EXRET_full(end - T + 1 : end , :);
EXRET_ll = EXRET_ll_full(end - T + 1 : end , :);

exret = exret_full(end - T + 1 : end , :);
exret_ll = exret_ll_full(end - T + 1 : end , :);

%% 4. STORE SIMULATIONS
SIMDATA = cell(length(Var_Names) , 1);
for i = 1 : length(Var_Names)
    SIMDATA{i} = eval(Var_Names{i});
end

end